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Abstract 

We present a new technique for a numerical analysis of the phase structure 
of the 2D Hubbard model as a function of the hole chemical potential. The 
grand canonical partition function for the model is obtained via Monte Carlo 
simulations. The dependence of the hole occupation number on the chemical 
potential and the temperature is evaluated. These calculations, together 
with a study of the Yang-Lee zeros of the grand canonical partition function, 
show evidence of a phase transition at zero temperature and particle density 
below half-filling. The binding energy of a pair of holes is calculated in the 
low temperature regime and the possibility for pairing is explored. 



1. Introduction 

For many years the Hubbard model, and other related systems with a finite 
density of electrons, has attracted much attention in the field of numeri- 
cal simulations The main interest in these simulations arose from 
the suggested relation between the planar Hubbard model and high Tc su- 
perconductivity (HTSC) 0, Since the model is non-relativistic, its 
analysis avoids some of the problems associated with relativistic fermions, 
such as fermion doubling p[. However, the inherent difficulties of simulating 
fermions at finite density remain: integration over the fermionic degrees of 
freedom leads to a non-positive integration measure in the path integrals. 
This arises from the determinant of the fermion matrix being non-positive 
definite and any importance sampling based on a partition function propor- 
tional to this determinant loses its normal meaning. 

In the Hubbard model, this problem manifests itself via the so-called sign 
problem of the partition function. The partition function of the half-filled 
Hubbard model (one electron per lattice site) is always positive, being a 
product of a determinant of a matrix and its hermitian conjugate. However, 
if one introduces a finite real chemical potential for the impurities (holes), 
then this positivity is lost, and configurations with real determinants, but 
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with negative signs can occur. Unfortunately, the finite density of holes is 
the case of greatest physical interest, as the impurities play an essential role 
in the superconducting transition 

In this paper we propose an analysis of the Hubbard model which treats 
the fermion dynamics in a rigorous manner. This method is similar to one 
applied to the chiral phase transition at finite density QCD, and is based 
on a study of the zeros of the grand canonical partition function in the 
complex fugacity (e'^^) plane [0, p. Here we simulate the configurations of 
the Hubbard- Stratonovich (HS) fields (given by Ising variables) at half-filling, 
as well as at finite doping fraction and expand the grand canonical partition 
function (GCPF) as a polynomial in e^'^ (or equivalently in e^^ + e~^'^)[0]. 
The coefficients of this polynomial are averaged over an equilibrated ensemble 
of configurations. The distribution of the zeros of this polynomial in the 
neighbourhood of the physical region, /x real, can indicate a phase transition. 
In particular, their scaling behaviour with respect to the lattice size can 
indicate if real zeros can occur in the infinite lattice limit. Such zeros would 
correspond to divergences within the theory at the corresponding values of 
the fugacity. 

The simulation is performed on a 2D spatial square lattice with the third 
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euclidean time dimension corresponding to the inverse temperature (3. The 
updating procedure is the one described by White et al. [||. 

As we show below, the evaluation of the partition function as an explicit 
polynomial in the fugacity variable permits analysis of the superconducting 
properties of the model for various values of the hole density. 

In Section 2 we describe the construction of the partition function of the 
Hubbard model as a polynomial in the fugacity variable. Section 3 sum- 
marizes the simulation procedure and the the measurements, and in Section 
4 we present the preliminary numerical results for the critical value of the 
chemical potential. 

2. Finite-temperature partition function 

The original Hubbard hamiltonian is given by: 

H = -tY, 4,aCj,^ + U Y.^ni+ - - ]-) - fiY,{ni+ + n^.) (1) 

i,j,a i i 

where the i,j denote the nearest neighbour spatial lattice sites, a is the 
spin degree of freedom and rii^ is the electron number operator cJ^Qo-. The 
constants t and U correspond to the hopping parameter and the on-site 
Coulomb repulsion respectively. The chemical potential /i is introduced such 
that fi = corresponds to half-filling. 
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The finite temperature grand canonical partition function (GCPF) is 
given by: 



wliere (3 is tlie inverse temperature. 

Tlie finite temperature is represented on the lattice by extending the 
spatial lattice in the imaginary-time direction and relating the inverse tem- 
perature (3 to the length of the time dimension rir hy (3 = rirdt, where dt 
is the length of the time step. Following Hirsch Jl| and White et al. 0] we 
rewrite the partition function as: 



where K corresponds to the nearest neighbour hopping term in the Hubbard 
hamiltonian (0) and V to the onsite interaction including quartic products 
of fermion fields. This decomposition, based on the Trotter formula [^], 
introduces a systematic error proportional to dt^. The quartic interaction 
can be rewritten in terms of Ising fields Si^i using the discrete Hubbard- 
Stratonovich transformation |]T|: 



(2) 



Z = Tr{e 



—dtV —dtK dtii\nr 



(3) 




dtU 1 



dtSi^l\{:ni+-ni-) 



(4) 
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where i, I is the space-time index of a lattice site and the couphng A is related 
to the original on-site repulsion constant by: 

cosh (dtX) = exp (~^)- (5) 

This linearization of the interaction enables one to integrate out the 
fermionic degrees of freedom and the resulting partition function is written 
as an ensemble average of a product of two determinants: 

Z = z = det{M+)det{M-) (6) 

such that 

M± = (/ + p±) = (/ + n Bh (7) 

1=1 

where the matrices Bf^ are defined as 

with Vij = SijSij and K the matrix connecting nearest-neighbours sites with 
the hopping parameter t = 1. The matrices in (|^) and (|^) are of size {nxHy) x 
[rixny), corresponding to the spatial size of the lattice. However, det{M^) 
can be represented as a determinant of an {rixnyn-r) square matrix of the 



6 



form 



( I 

B 








/ / 



(9) 



V ... -Bl 

a fact which is exploited below in the evaluation of the partition function Z. 

The expectation value of a physical observable at chemical potential /i, 
< O is given by: 



<0>, 



JOzifi) 



(10) 



where the sum over the configurations of Ising fields is denoted by an integral 
Since z{fi) is not positive definite for Re(/i) 7^ we weight the ensemble of 
configurations by the absolute value of at some = fiQ. Thus 



<0>, 



>, 



The partition function Z{^) is given by 

f(/^o)| 



(11) 



MO 



(12) 



The normalization of the GCPF is irrelevant as can be seen from eg . (|lTl) . 



The particle-hole transformation ||1|,||1CI|| 



din 



(13) 
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is equivalent to the unitary transformation 



/i ... 0\ /-i ... 

-i ... + i .. 

i ... ' -i .. 

\0 ... -i) \ ... i J 



(14) 



which reverses the sign of the hopping term K. Hence on an even sized lattice 
the determinant of e'^*^ is 1. Applying this transformation to the statistical 
weight gives: 



Equation (|15]) suggests two different ways of expansion of the partition 

function. The first one based on 

nA,((e''^ + e"^^) + (A, + ^)) X e'""-"^^ 



and the second on: 



n=0 



(16) 



(17) 



where the Aj are the eigenvalues of the matrix P\^_^- Note that the expansion 
coefficients 6„ are the canonical partition functions of the n-electron excita- 
tions above half-filling, (with n < corresponding to holes). We note here 



that eqns.(0),(|l7|) follow from the fact that the eigenvalues of P~_^ are either 
real or appear in complex conjugate pairs. The coefficients of the character- 
istic polynomials, namely {a„} and {bn}, are obtained from ([T^) , (P^) by the 
recursion procedure described in The sign problem manifests itself in 



the fluctuating signs of these coefficients from configuration to configuration 
of equilibrated Ising fields. The expansion coefficients for a grand canonical 
partition function (GCPF) are then obtained by averaging the coefficients 
of each of these polynomials over the ensemble of configurations. A similar 
procedure has been applied in the study of the chiral phase transition in fi- 
nite density lattice QCD and in the evaluation of the critical mass in lattice 
QCD 0. 

At /io = (and at any imaginary chemical potential 0]) z{fio) is clearly 
positive. With alternative choice of the updating fiQ the GCPF, Z{^o), is 



equal to the average sign of the weight function z(/io)[|T2[,|[T^. We have 
performed calculations using updating at half-filling and at /xq 0. The 
latter choice of the updating chemical potential is important for simulations 
performed at low temperatures. We will show below that it provides re- 
sults consistent with the half-filling updating at relatively high temperatures 
(/3 < 5.) while at higher values of (3 it provides better numerical stability in 
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obtaining the expansion coefficients of the GCPF, corresponding to the finite 
hole occupation. 

As the temperature is lowered the bounds on the eigenvalues of the matrix 
P|^_^, which are found via the Lanczos algorithm, become very large. Simu- 
lating configurations at /3 = 10., dt = 0.125 we need to handle a lattice with 
Ut = 80. For this set of parameters we find that the eigenvalues vary between 
10^^ and 10~^^ on a 10^ lattice which damages severely the efficiency and the 
accuracy of the whole calculation. However, the characteristic polynomial 
can be also obtained from the determinant of the the {nrUxny) x {nrUxny) 
matrix / + given in (|^), using the eigenvalues of 

It follows from the structure of A~^_^ that its eigenvalues have a sym- 
metry: if Aj is an eigenvalue so is A^e [n = 1, ...rir — 1)- The coefficients 
of the z expansion in the fugacity powers are actually functions of the A"^. 
The variation in these eigenvalues is significantly smaller, but the matrix 
to be diagonalized is n"^ times bigger, leading to a more time consuming 
diagonalization procedure. 

The method used in our calculations consists of representing the deter- 
minant of d^) as 

det{I + A) (18) 
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with 



/ ... BiB2...Bn,\ 

-Bnt+lBnt+2---- ... 



A 



(19) 



V -...5„_iB„^ / 

where the matrix A is of the size {nxUy'^) x (uxny^). (In the last equation 
the B matrices are taken at /i = 0.) The eigenvalues of this matrix have a 
reduced symmetry Z^^ -^n^/nt and thus variations of a larger magnitude 
than those of A, but the diagonalization procedure is more efficient. On the 
other hand its eigenvalues are varying in a smaller range than the eigenvalues 
of the total time ordered product nr=i ^i- By finding the most appropriate 
ratio riT-fnt we succeed to obtain the eigenvalues of A with the required 
precision and then taking the n^/rit power of them we obtain the {Aj} and 
evaluate the expansions (|1|) and (0). Note that since the coefficients of 



the characteristic polynomials for ([T6[) and ([17|) depend only on Aj + A~^, the 
above procedure, although introducing large variations in the Aj's themselves, 
does not lead to significant errors in the coefficients around half-filling. It is 
these coefficients which determine the behaviour of the smallest zeros of the 
GCPF and thus the phase structure of the model. 

The physical observables are derived using eq.(|lT]). For a given config- 
uration of Ising fields {sj ^}, the value of an operator O can be calculated 
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as a polynomial in the fugacity variable. Knowing the coefficient of each 
power of the fugacity in the z expansion, one can then easily construct the 
corresponding coefficients for the contribution of this term to the observable 
by averaging each coefficient over the equilibrated ensemble. In this paper 
we measure only the expansion of the GCPF and hence we can predict the 
critical value of the chemical potential for which the relative weight of the the 
state with a finite doping fraction will be of the same order as the half-ffiled 
state. We expect that this prediction will be reflected in the behaviour of the 
hole occupation number and consistent with the critical values of /i obtained 
from an analysis of the complex zeros of the GCPF. Moreover, following the 
suggestion of Dagotto et al.0 we evaluate the energy gap between the one 
pair state and the half-filled state. 

3. Simulations and measurements 

The simulation procedure is based on the algorithm described by White et 
al. Q . We use Metropolis algorithm for updating the configurations of Ising 
variables. Here we describe the procedure for the half-filling updating. (For 
clarity, we omit the spin labels in the following.) The generalization for the 
finite /xq updating is straightforward. 

The simulation starts from an arbitrary configuration for which we cal- 

12 



culate the equal-time Green's function on a time slice / 

G{1) = {I + P{1))-' (20) 

where P{1) is a time ordered product of the form 

Bi...BiBn^...Bi^i (21) 

To reduce the numerical errors in the evaluation of these products we apply 
the modified Gram-Schmidt (MGS) decomposition as proposed by White 
et al. [Q. Decomposing the products of each four matrices in (^Tj) into a 
product of an orthogonal matrix, a diagonal one and an upper-triangular 
matrix whose diagonal elements equal to one enables us to deal with large 
variations in the matrix elements. The inversion of the (J + P(/)) appearing 



in the r.h.s. of ( |2(]|) is also simplified by the MGS procedure. Once the equal- 
time propagator is evaluated on a time-slice / we flip one of the spins Si^i and 
and accept the new configuration with respect to the ratio of the new and 
old statistical weights, defined as: 

detM'^ . 

where 5{si^i) is the difference in the potential term due to the flipped spin. 

This ratio is determined by the value of the equal-time Green's function 
(by its diagonal term corresponding to the flipped spin) and by a matrix A 
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with only one nonzero element: 

A(^,0,,fc = e"2'^*"^-'5„-5,,fc (23) 

If the new configuration is accepted we calculate the new Green's function 
G{iy corresponding to this configuration using: 

G{iy = G{1) - 4 X G{l)A{i, - G{1)) (24) 

with 

detM' 

R=^^ = l + {l-G,,)/\{t,l)u (25) 

The nonlocal impact of the updated configuration is represented in the new 
Green's function by eq. (p^). 

When the updating of the Ising fields on the 1-th time-slice is completed 
we move to the next time slice using the relation: 

G{1 + 1) = fi,+iG'(/)fir+i- (26) 

Following the suggestion of |Q, we evaluated the Green's function from 
scratch every four time steps. This procedure is very time consuming, but 
is required in order that the numerical errors accumulated using eq.(p6D are 
kept under control. Since the hopping term in the Hamiltonian is constant, 
the construction of the Green's function from scratch reduces to a redefini- 
tion of the interaction matrix e'^*^ on each time-slice due to the updated Ising 
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fields. The computational effort involved in this procedure is minor as is a 
diagonal matrix and its exponentiation is fast. The hopping part is exponen- 
tiated only once at the start of the simulation procedure. Instead of using the 
checkerboard decomposition suggested by White et al. we expanded e"*^*^ 
taking advantage of its sparseness. Since this expansion is performed only 
once it can be done up to an arbitrary high order. We checked that taking 
a 10-th order expansion provided us with sufficiently high precision for the 
parameters used in the simulations described below. 

4. Results and conclusions 

The expansion coefficients of the GCPF as a polynomial in the fugacity vari- 
able were calculated at several values of the inverse temperature (3 with the 
Coulomb repulsion fixed aX U = At (see eq.(|l|). We present results obtained 
from simulations performed at half-filling and at chemical potential /io = 0.9. 
The spatial size of the lattice is 4^ throughout apart from one simulation at 
half-filling on a 10^ lattice at /3 = 10.0. 

Half-filling results: The half-filling simulation was performed at /3 = 
0.3,1.2,2.5,5.0 and 10.0 with = 4,16,20,40 and 80 respectively. The 
euclidean time-spacing dt was varied to check the numerical stability of the 



simulations and to allow comparison with the results of other groups [|T4 
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The number of configurations required to get sufficiently low errors in the 
expansion coefficients is, in general, greater than 2000 and increases with /3. 
The particle-hole symmetry is manifested via the equality between the coef- 
ficients of e'^"'^ and e"'^"'^. This symmetry follows directly from eqs. (|l6| , p!7|) . 
The Zn =< bn > coefficient corresponds to the contribution of the n-hole 
state to the GCPF (canonical parition functions) [§] and Zq is the canonical 
partition function for the half-filled state. The latter is obtained with low 
error after a small number of measurements. The higher order coefficients 
require averaging over a larger number of configurations. For /? < 5 all the 
2{nxnyY + 1 averaged coefficients were found to be positive. 

As we extend our treatment to larger values of /3, negative coefficients 
appear in the expansion of the GCPF, but with large errors. Note that the 
coefficients arising from a single configuration do not have to be positive, as 
only the ensemble averages are identified as the canonical partition function 
at a given particle number. However, the low temperature simulation at 
half-filling does lead to a high variation of the coefficients in different field 
configurations and thus to the high errors. This is because updating at 
/i = minimizes the fluctuations in only zq (half-filling) which dominates 
the statistical weight at high (3. The large variations in these coefficients 
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indicates that this statistical weight, namely the determinant at /io = 0, 
becomes inefficient as /3 gets large. A more appropriate choice is to update 
at jj,Q 7^ 0. 

Simulations at /iq 7^ 0: We performed our simulations at /io = 0.9 
and updated with respect to the absolute value of the weight function as 
described in the previous section at /5 = 2.5, 3., 5.0, 5.4, 6.0 and 7.5 with rir = 
20, 40, 40, 72, 48 and 60 respectively. In this simulation the value of the GCPF 
at /io is the average sign of det{M~^)det{M~) . This requirement provides a 
useful check as to our numerical accuracy in extracting the coefficients of the 
characteristic polynomial. At (3 = 2.5 and = 5. we compared the results 
of these simulations with those performed at half-filling updating and found 
that the normalized coefficients are equal within the statistical error. 

In Table 1 and Fig.l we present the expansion coefficients based on eq. 
( P^Tf ) obtained from simulations performed either at half-filling or at updating 
chemical potential /io = 0.9. The normalization of the GCPF is chosen such 
that Z(/i = 0) = 1. Our results are consistent with those obtained by Moreo 
et a/.0. 

Fig.l shows the coefficients Zn {n = 0,1,2,3) as a function of (3, with 
normalization such that the GCPF Z(/i = 0) = 1. The general tendency is 
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a sharp decrease of the coefficients with higher occupation number. 

Figs. 2 and 3 show the hole density and the susceptibihty respectively, as 
a function of the chemical potential. In these exploratory simulations, the 
coefficients corresponding to high occupation number are determined with 
large errors. However, the peak in the susceptibility at /i ~ 1.0 is due to the 
low occupation levels which are determined with small errors. The structures 
at fi > 1.2 do depend on the levels with large errors and require further 
investigation. As the temperature decreases the peak in the susceptibility 
sharpens significantly in the region 0.75 < fi < 1.25, especially for f3 > 5.0. 
At that P the susceptibility does seem to signal some change in behaviour. 
This could be associated with the onset of a phase transition related to the 
occupation of holes. We explore this possibility further by analyzing the 
zeros of the GCPF in the complex fugacity plane. 

We do this by finding the zeros of the averaged polynomial eq.(^). Since 
the large n coefficients are evaluated with relatively low precision we checked 
the stability of the small zeros under truncation of the polynomial to ri = 
4. Table 2 gives the two smallest zeros for various values of the inverse 
temperature. In Fig.4 we plot these zeros of the partition function in the 
first quadrant of the complex ^ plane. The zeros have the symmetries that 
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if /X is a zero, so also is — and their complex conjugates. 

To a very good approximation, the imaginary part of these zeros of the 
partition function scales as ^. For any finite value of (3 the fugacity would 
remain negative yielding no phase transition in the physical region. However, 
at zero temperature {(3 — oo) the Im(/Xc) vanishes and a phase transition may 
be possible. To extrapolate to the low temperature behaviour of the zeros we 
note an almost linear dependence between the imaginary part of the lowest 
zero and its real part, for values of /3 > 5. The linear fit of these zeros is shown 
in Fig.5. To clarify this point we show in Fig.6 the hnear fit of Re(//) x /3 vs. 
P in the same region of (3. The lowest zero at /3 = 2.5 and 3 does not exhibit 
this scaling behaviour. Measuring the lowest zero of the partition function 
at high temperatures (/3 = 0.3 and (3 = 1.2) we find that the real parts of the 
zero to be 3.0 and 0.721 respectively with imaginary part |. The above is 
consistent with the lowest zero, /ic, scaling such that Re(/Xc) x /? is constant 
in the high temperature regime. At lower temperatures, /? = 2.5 — 3, there 
is a crossover region into the low temperature regime where 

i?e/Xc = -2.5//3 + 1.1. (27) 

This qualitative distinction between the high and low temperature regimes 
arises from a clear difference in the relative contributions of the finite particle 
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number states to the grand canonical partition function (see Table 1). In the 
high temperature regime the first two canonical partition functions are of 
the same order as the half-filled level and hence the corresponding states are 
excited even at zero chemical potential. On the other hand, these states in 
the low temperature GCPF are only excited by non-zero chemical potential. 
Based on the above observations, we conclude that there is the possibihty of 
a phase transition at zero temperature and /ic ~ 1-1 but that no signal has 
been found for a phase transition at T > 0. Of course, the above simulations 
have been performed on a small system. There may well be large finite 
volume effects. 

The above conclusion - that there may be a phase transition at zero 
temperature - will not alter if the lowest zero in the fugacity plane remains 
real (and hence necessarily negative). If there is a phase transition at nonzero 
temperature, then some complex zeros must be in the complex fugacity plane 
with Re{e^^) > and pinch the positive real axis in the infinite spatial volume 
limit. No signal of the possible onset of this mechanism was observed on the 
4^ lattice, i.e. no zeros were found in the first or fourth quadrants of the 
complex e^^ plane. 

As we increase the spatial lattice size we should therefore observe either an 
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increase in the density of zeros adjacent to, or on, the negative real fugacity 
axis or, if the alternative mechanism is masked by finite size effects on the 
4^ lattice, a migration of zeros to the Re{e^^) > complex half-plane. It 
is also important to confirm the scaling behaviour of eq.(^) for the lowest 
zero. Zeros adjacent to this one should also scale in an analogous manner 
with the temperature so that there is a well defined locus of zeros in the zero 
temperature limit. We intend to extend our simulations to 6^ and 8^ spatial 
lattices. 

The absence of a critical positive fugacity above zero temperature can 



be interpreted in part as a realization of the Mermin- Wagner theorem [|I5 
which claims absence of magnetic ordering in two dimensional spin systems 
at non-zero temperature. This theorem is in particular relevant to HTSC 
models based on the Heisenberg antiferromagnetic Hamiltonian resulting in 
the strong coupling treatment of the Hubbard model. If the isotropy of 
Heisenberg antiferromagnet is violated, e.g. by an interlayer interaction, the 
conditions of the theorem do not hold, thus relating the vanishing critical 
temperature for the superconductivity to the isotropy of the effective nearest 



neighbours coupling[16|. The 2D Hubbard model with effective interlayer in- 



teraction was recently studied by M. Frick et al. using the Projector Monte 
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Carlo technique ||1 7|| , p8[] . Their results show some evidence for HTSC. We 
note that the analysis described in our work can also be applied to the ex- 
tended Hubbard model including effective interlayer interactions. 

Finally, following the suggestion of Dagotto et al. we calculate the 
binding energy of holes. The energy gap between the two and one-hole ground 
states, E2 — El, is given by the slope of the linear fit to the ratio log — vs. P 
(see Figs. 7, 8). Fig. 7 presents the fit for the data in the range 2.5 < (3 < 7.5 
while, in Fig. 8, we fit only the results for (3 > 5.0 in the light of the discussion 
above. The first fit shows an energy gap of 0.87 ± 0.02 which is very close to 
the result of Dagotto et al. (0.88 ± .0.02) [Q. The energy gap obtained from 
the second fit is 0.85 ± 0.02. Analogous fits, (Figs. 9, 10), for log — vs. (3 give 
the energy gap between the ground state with one hole and the corresponding 
state at half-filling, Ei-Eq. The fits give ^i-^o = 0.83±0.02 and 0.95±0.02 
respectively. The binding energy of holes is given by the difference 

(E2-^i)-(Ei-Eo). (28) 

Thus the first fit yields positive binding energy with no pair creation expected 
while the latter suggsts a binding energy of —0.1. Since this derivation of 
the energy gaps is valid only in the low temperature regime, which is distinct 
from the high temperature one, the second fit seems to be more appropriate. 
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We note that we have not included the spin wave contribution to the ratio of 
the canonical partition functions Z\ and -ZoO^ since at low temperature, the 
spin wave contribution should become small. The result of our low temper- 
ature fit is close to that obtained by Dagotto et al. {Ei — Eq = 0.98 ± 0.02). 
However, their fit included lower (3 data thus requiring a spin wave contribu- 
tion. Inclusion of a spin wave contribution to our canonical partition function 
at half-filling would raise the estimate of the one-hole ground state energy 
even higher, thus increasing the binding energy. 

We conclude with suggested extensions of the above method. The pre- 
dicted zero temperature phase transition can be confirmed by a lower tem- 
perature study. With this in mind, simulations at = 12 (jIt = 160) are 
under current investigation. A study of the finite size effects, in particular, 
the scaling properties of Im(/ic) as a function of the volume is also necessary. 
It is also important to perform longer simulations in order that the analysis 
can be extended to larger values of the chemical potential. One can also gen- 
eralize the method described above to derive polynomial expansions in the 
fugacity variable for other physical observables and thus extend the study of 
the nature of the phase transition and of its possible relevance to high Tc 
superconductivity. This study would involve examining the persistence of 
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the antiferromagnetic order into the finite doping region p9|]. 
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Table 1. 





Zq 


Zl 


Z2 


z?. 




coeff. 


err. 


coeff. 


err. 


coeff. 


err. 


coeff. 


err. 


0.3 


0.171 


2.E-4 


0.156 


1.5E-4 


0.119 


3.E-4 


0.074 


6.E-5 


1.2 


0.297 


0.001 


0.223 


2.E-4 


0.097 


4.E-4 


0.025 


2.5E-4 


2.5 


0.477 


0.002 


0.220 


0.001 


0.039 


5.E-4 


0.003 


8.E-5 


3. 


0.575 


0.016 


0.189 


0.002 


0.022 


l.E-4 


0.001 


l.E-5 


5. 


0.872 


0.005 


0.063 


0.002 


0.001 


2.E-4 


5.E-6 


3.E-6 


5.4 


0.901 


0.061 


0.049 


0.002 


7.E-4 


3.E-5 


4.E-6 


2.E-7 


6. 


0.946 


0.056 


0.027 


0.001 


2.E-3 


l.E-5 


6.E-7 


4.E-8 


7.5 


0.986 


0.063 


0.007 


5.E-4 


2.E-5 


l.E-6 


l.E-8 


l.E-9 
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Table 2. 





Re/i 


Im/i 


(Im/i) X P 




Full 


Trunc. 


Full 


Trunc. 


Full 


Trunc. 


2.5 


0.526 


0.464 


1.257 


1.257 


3.141 


3.141 




0.821 


0.731 


1.257 


0.952 


3.141 


2.380 


3. 


0.582 


0538 


1.047 


1.047 


3.141 


3.141 




0.764 


0.749 


1.047 


0.819 


3.141 


2.457 


5. 


0.609 


0.607 


0.628 


0.628 


3.140 


3.140 




0.775 


0.818 


0.628 


0.628 


3.140 


3.140 


5.4 


0.627 


0.625 


0.582 


0.582 


3.143 


3.143 




0.884 


0.820 


0.572 


0.477 


3.089 


2.576 


6. 


0.685 


0.680 


0.524 


0.524 


3.144 


3.144 




0.801 


0.848 


0.524 


0.453 


3.144 


2.718 


7.5 


0.768 


0.785 


0.380 


0.397 


2.850 


2.978 




0.990 


0.872 


0.310 


0.419 


2.325 


3.142 
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Figure Captions 

1. Fig.l Zo,Zi,Z2, and as a function of Their associated statistical 
errors are also shown if larger than the resolution. 

2. Fig. 2 Average hole density (doping fraction) as a function of n for P 
between 2.5 and 7.5: 

(3 = 2.5 and 3 — solid lines 

P — 5.0 — dashed line 

P = 5.4 — dashed-dotted line 

P — 6.0 — dotted hne 

P = 7.5 — solid line; 

3. Fig. 3 Susceptibility as a function of /i for P between 2.5 and 7.5: 
P — 2.5 and 3 — solid lines 

P — 5.0 — dashed line 

P = 5.4 — dashed-dotted line 

P — 6.0 — dotted hne 

P = 7.5 — solid line; 

4. Fig. 4 Two smallest zeros of the GCPF in the first quadrant of the 
complex II plane for P between 2.5 and 7.5: 
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/? = 


2.5 


— white dotted squares 




3.0 


— black crosses 


/? = 


5.0 


— white circles 




5.4 


— black circles 


= 


6.0 


— black squares 


/3 = 


7.5 


— white squares; 



5. Fig. 5 Im(/i,c) plotted against Re(/Xc)- The solid line is a linear fit for /5 
between 5.0 and 7.5; 

6. Fig. 6 Re(/ic)/3 as a function of /3. The solid line is a hnear fit for beta 
between 5.0 and 7.5; 

7. Fig. 7 log(^2/^i) as a function of The solid line is a linear fit for 
2.5 <f3< 7.5; 

8. Fig. 8 log(2;2/2;i) as a function of (3. The solid line is a linear fit for 
5.0 < /3 < 7.5; 

9. Fig. 9 log(2;i/2;o) as a function of /3. The sohd fine is a linear fit for 
2.5 < (3 < 7.5; 

10. Fig. 10 log(2;i/2;o) as a function of The solid line is a linear fit for 
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5.0 < /3 < 7.5; 
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Table Captions 

1. Table 1 Expansion coefficients of the GCPF for (3 between 0.3 and 2.5. 
zq corresponds to the half-filling and z„ to the n— hole coefficients. 

2. Table 2 Smallest zeros of the GCPF in the complex ji plane for (3 
between 2.5 and 7.5 . The zeros are obtained from the full polynomial 
and the truncated one (up to 4-hole coefficient), to check the numerical 
stabihty of the results. 
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